program drop _all
program define PANELpbias, rclass
version 13.1

syntax, studies(integer) estperstudy(integer) totalobs(integer) alpha(real) ///
theta(real) obs(integer)
        
// Remove existing variables
drop _all

//We first create the matrix to store the results of each study
set matsize 5000
matrix A = J(`totalobs',1,.)
matrix B = J(`totalobs',1,.)
matrix C = J(`totalobs',1,.)
matrix D = J(`totalobs',1,.)
		
forvalues i = 1/`studies' {
scalar lambdai = 0.5+30*runiform()
scalar alphai = `alpha'+2*rnormal()
forvalues j = 1/`estperstudy' {
// STEP ONE: Create the data for each study and estimate an effect
clear
set obs 100
generate x = rnormal()
// Note that each "study" has a difference error variance, causing the estimate 
// of the effect to be estimated with varying degrees of precision
// This is for random effects
scalar lambdaij = lambdai+`theta'*runiform()
scalar alphaij = alphai+0.5*rnormal()
generate e = lambdaij*rnormal()
generate y = 1 + alphaij*x + e
// This is for fixed effects
// scalar lambda = 0.2+30*runiform()
// generate e = lambda*rnormal()
// generate y = 1 + `alpha'*x + e
quietly regress y x 
scalar coef = _b[x]
scalar secoef = _se[x]
scalar tcoef = coef/secoef
scalar ID = `i'
scalar obsno = (`i'-1)*`estperstudy'+`j'
	
// First run this program once to get the pre-publication study sample data
// To get post-publication study sample data, uncomment one of the two sections
// below.
matrix A[obsno,1] = coef
matrix B[obsno,1] = secoef


matrix C[obsno,1] = tcoef
matrix D[obsno,1] = ID
}
}

// The next set of commands moves the data out of matrices and reformats them as 
// standard Stata data series.  We have now completed generating our individual
// studies and we now move into the "meta-analysis" stage.
matrix bob = A,B,C,D
svmat bob
rename bob1 effect 
rename bob2 seeffect 
rename bob3 teffect 
rename bob4 ID
generate pet = (1/seeffect)

// This set of commands imposese the publication bias, where the publication criterion
// is either that (i) the t-stat must be greater than or equal to 2, or (ii) the 
// estimated effect is positive.  The commands below implement the assumption that a
// study must have at least 7 out of 10 estimates that satisfy the publication 
// criterion in order for the study to be "published."
generate dummy = 1
//replace dummy = 0 if abs(teffect) < 2
replace dummy = 0 if effect < 0
by ID, sort: egen select = mean(dummy)
// Not sure why this happens, but if I put select<0.70, it kicks out the studies that
// have 7 estimates that satisfy the publication criterion.  So I set select<0.65.  
// Studies are omitted from the "meta-analysis" sample by replacing the relevant variables
// with missing values.
replace effect = cond(select<0.65,.,effect)
replace seeffect = cond(select<0.65,.,seeffect)
replace teffect = cond(select<0.65,.,teffect)
replace pet = cond(select<0.65,.,pet)
gen fepet = pet
gen feteffect = teffect
	
// This creates dummy variables for each of the 100 studies
// The dummy variables take names dum1 to dum100
tab ID, gen(dum)
	
// This creates study-specific SE terms for use in the PEESE
// according to equation 5.7 on page 85 of S&D
forvalues i = 1/100 {
generate SE`i' = seeffect*dum`i'
}
	
		
// This estimate produces the PET version of the "publication bias"- 
// corrected effect estimate
// NOTE1:  According to equation 5.6 on page 85 of S&D, the bias-corrected
// effect is given by the coefficients on the respective precision terms, pet*.
// The specification below forces all the effects to be the same, while
// allowing for fixed effects to correct for bias-associated with estimate SEs.
// NOTE2: Also note that while all the dummy variables will not be included in the
// meta-analysis sample, this is not a problem because STATA will automatically
// kick out the ones that don't belong.
regress teffect dum1-dum100 pet, vce(cluster ID)
return scalar effect_PET = _b[pet]
scalar effect_PET = _b[pet]
test _b[_cons] = 0
return scalar pvalue_FAT = r(p)
scalar pvalue_FAT = r(p)
test pet = `alpha'
return scalar pvalue_PET = r(p)
scalar pvalue_PET = r(p)
test pet = 0
return scalar pvalue_PETFPP = r(p)
scalar pvalue_PETFPP = r(p)
		
// This estimate produces the PEESE version of the "publication bias"- 
// corrected effect estimate.  It is based on equation 5.7 on page 85 of S&D.
// See notes from above.
regress teffect SE1-SE100 pet, noc vce(cluster ID)
return scalar effect_PEESE = _b[pet]
scalar effect_PEESE = _b[pet]
test pet = `alpha'
return scalar pvalue_PEESE = r(p)
scalar pvalue_PEESE = r(p)
	
	
// This estimate produces the WLS-FE estimate of the effect
regress feteffect fepet, noc vce(cluster ID)
return scalar effect_FE = _b[fepet]
test fepet = `alpha'
return scalar pvalue_FE = r(p)	
	
		
// This estimate produces the WLS-RE estimate of the effect
// NOTE: We use the Method of Moments (mm) option of metareg
// because the maximum likelihood procedure had too many instances
// of failure to optimize.  Method of Moments does not require
// iteration and thus avoids this problem.
quietly metareg effect, wsse(seeffect) mm
scalar tau2 = e(tau2)
gen revarR= seeffect^2 + tau2
gen reseR = sqrt(revarR)
gen reteffect = effect/reseR
gen repet = 1/reseR
regress reteffect repet, noc vce(cluster ID)
return scalar effect_RE = _b[repet]
test repet = `alpha'
return scalar pvalue_RE = r(p)	
	
return scalar effect_FPP = effect_PET
return scalar pvalue_FPP = pvalue_PET


if pvalue_PETFPP < 0.05 {
return scalar effect_FPP = effect_PEESE
return scalar pvalue_FPP = pvalue_PEESE
}	



end
